//
// RandomSpin.cxx
//
// Created by TB on 11/23/11.
// Copyright 2011 BNL. All rights reserved.
//

#include <map>

#include <TRandom3.h>

#include "tppmc/Enumerations.h" // For LHAPDF::Flavour
#include "tppmc/Particle.h"
#include "tppmc/PartonContent.h"
#include "tppmc/PartonDistribution.h"
#include "tppmc/RandomSpin.h"

namespace {
   
   // Convert PDG code for a particle to the corresponding LHAPDF::Flavour.
   // Returns -999 if the input PDG code is not a parton/gluon/photon.
   // Turn this utility function into a class that other files
   // can access, as it would be handy.
   int pdgToLhapdf(int pdgCode) {
      
      static std::map<int, int> mTable;
      
      if(mTable.empty()) {
         int pdg[14] = {
            1,  2,  3,
            4,  5,  6,
            -1, -2, -3,
            -4, -5, -6,
            21, 22
         };
         int lhapdf[14] = {
            LHAPDF::DOWN,  LHAPDF::UP,     LHAPDF::STRANGE,
            LHAPDF::CHARM, LHAPDF::BOTTOM, LHAPDF::TOP,
            LHAPDF::DBAR,  LHAPDF::UBAR,   LHAPDF::SBAR,
            LHAPDF::CBAR,  LHAPDF::BBAR,   LHAPDF::TBAR,
            LHAPDF::GLUON, LHAPDF::PHOTON
         };
         for(unsigned i(0); i < 14; ++i) {
            mTable.insert(std::make_pair(pdg[i], lhapdf[i]));
         } // for
      } // if
      
      // Find the entry
      const std::map<int, int>::const_iterator iter =
      mTable.find(pdgCode);
      
      // Didn't find it:
      if(mTable.end() == iter) return -999;
      
      return iter->second;
   }
}

namespace tppmc {
   
   
   RandomSpin::RandomSpin(PartonDistribution* unpolarised,
                          PartonDistribution* polarised)
   : mUnpolarised(unpolarised)
   , mPolarised(polarised) {
   }
      
   
   RandomSpin::~RandomSpin() {
   }
   
   int RandomSpin::Generate(const pythpart& parton,
                            double scale) const {
      // Evalute the parton distributions for this kinematic point
      PartonContent unpol;
      PartonContent pol;
      mUnpolarised->Evaluate(parton, scale, unpol);
      mPolarised->Evaluate(parton, scale, pol);
      
      const int id = pdgToLhapdf(parton.Type());
      if(-999 == id) return 0.;
      
      // Get the PDF values for this parton type:
//      const double unpol_pdf(unpol[id]);
//      const double pol_pdf(pol[id]);
      
      double prob_spin_up(0.);
      
      if(unpol[id] not_eq 0.) {
         prob_spin_up = 0.5 * (1. + (pol[id] / unpol[id])); 
      } // if
      else {
         prob_spin_up = 0.5;
//         std::cerr <<
//         "WARNING: unpol_pdf=0.0, assuming prob_spin_up = 0.5. x = "<<
//         parton.X() << " Q2 = " << scale << " type = " << parton.Type() <<
//         std::endl;
      } // else
      
      // Determine the probability spin is positive
      // (along hadron spin) or opposite. 
      
      if(gRandom->Uniform(0., 1.) < prob_spin_up) {
         return 1; 
      } // if
      else {
         return -1;
      } // else
   }
   
} // namespace tppmc
